function ComputeFiso()

dwi_pathname = 'C:\Ziang\Dataset\HumanConnectome\792564_processed\1X_Filtered\DWI_Filtered.raw';
gradient_pathname = 'C:\Ziang\Dataset\HumanConnectome\792564_processed\1X_Filtered\gradient_B3000.txt';
mask_pathname = 'C:\Ziang\Dataset\HumanConnectome\792564_processed\1X_Filtered\Brain_MASK.raw';
result_pathname = 'C:\Ziang\Dataset\HumanConnectome\792564_processed\1X_Filtered\Fiso.raw';

width = 114;
height = 150;
depth = 110;
n = 108;
n0 = 18;
b_value = 3000;

[DWI, orientation] = ReadDWI(dwi_pathname, gradient_pathname, n, width, height, depth);
MASK = ReadMask(mask_pathname, width, height, depth);

% compute Fiso for each voxel
Fiso = zeros(width, height, depth);
for z = 1 : depth
    fprintf('Slice z = %d ...\n', z);
    for y = 1 : height
        for x = 1 : width
            if MASK(x, y, z) > 0.5
                S =  DWI(n0+1:n, x, y, z);
                S0 = mean(DWI(1:n0, x, y, z));
                if (max(S) > S0)
                    Fiso(x, y, z) = 1.0;
                else                   
                    temp = S / S0;
                    d = min(log(temp));
                    d = -d / b_value;
                    Sa = mean(S);
                    Smax = S0 * sqrt(pi / (4 * b_value * d)) * erf(sqrt(b_value * d));
                    Smin = S0 * exp(-b_value * d);
                    Fiso(x, y, z) = (Smax - Sa) / (Smax - Smin);
                end
            else
                Fiso(x, y, z) = 1.0;
            end
            
        end
    end
end

% save result
SaveResult(result_pathname, Fiso, width, height, depth);
end

function SaveResult(pathname, Fiso, w, h, d)

fid = fopen(pathname, 'w');
for z = 1 : d
    for y = 1 : h
        fwrite(fid, Fiso(:, y, z), 'float');
    end
end
fclose(fid);


end
